Inaccuracy tests

Basically, how bad of an approximation can we do to the interpolation and still get away with 10% accuracy to the grid points? 10K, 0.05 dex in logg and [Fe/H] ?

As a comparison we need a high quality spectrum generated using our normal interpolation methods.


In [1]:
%matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as sp
from IPython.html import widgets # Widget definitions
from IPython.display import display, clear_output, HTML # Used to display widgets in the notebook
from IPython.html.widgets import interact, interactive, fixed


Using matplotlib backend: Qt4Agg

In [2]:
from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import TRES, HDF5Interface

myDataSpectrum = DataSpectrum.open("../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
myInstrument = TRES()

#myHDF5Interface = HDF5Interface("../libraries/PHOENIX_objgrid6000.hdf5")
myHDF5Interface = HDF5Interface("../libraries/PHOENIX_objgrid4000.hdf5")

myModel = Model(myDataSpectrum, myInstrument, myHDF5Interface, stellar_tuple=("temp", "logg", "Z", "vsini", "vz", "logOmega"), 
                cheb_tuple=("c1", "c2", "c3"), cov_tuple=("sigAmp", "logAmp", "l"), region_tuple=("h", "loga", "mu", "sigma"))
myOrderModel = myModel.OrderModels[0]


ss is len 2755
Creating OrderModel 0

In [5]:
from StellarSpectra.model import ModelHighAccuracy

#myHDF5InterfaceHighAccuracy = HDF5Interface("../libraries/PHOENIX_submaster.hdf5")
myHDF5InterfaceHighAccuracy = HDF5Interface("../libraries/PHOENIX_LkCa15.hdf5")
#myHDF5Interface = HDF5Interface("../libraries/PHOENIX_M.hdf5")

myModelHA = ModelHighAccuracy(myDataSpectrum, myInstrument, myHDF5InterfaceHighAccuracy, stellar_tuple=("temp", "logg", "Z", "vsini", "vz", "logOmega"), 
                cheb_tuple=("c1", "c2", "c3"), cov_tuple=("sigAmp", "logAmp", "l"), region_tuple=("h", "loga", "mu", "sigma"))
myOrderModelHA = myModelHA.OrderModels[0]


Grid stretches from 5134.978696599707 to 5235.9293774845055
wl_FFT is 0.0480588375883031 km/s
Creating OrderModel 0

In [3]:
myHDF5Interface.bounds


Out[3]:
{'alpha': (0.0, 0.0),
 'temp': (5900, 6100),
 'logg': (4.0, 5.0),
 'Z': (-0.5, 0.0)}

In [6]:
spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]

In [7]:
def plot_diff_spec(wl, spec0, spec1, spec2=None, norm=True, abs_err=False):
    '''
    Plot the difference between two spectra, if normalize, divide each by their mean.
    '''
    
    if norm:
        spec0 = spec0/np.mean(spec0)
        spec1 = spec1/np.mean(spec1)
        if np.all(spec2) is not None:
            spec2 = spec2/np.mean(spec2)
        
    fig, ax = plt.subplots(nrows=2, figsize=(11,8), sharex=True)
    ax[0].plot(wl, spec0, "b", label="spec0")
    ax[0].plot(wl, spec1, "r", label="spec1")
    if np.any(spec2):
        ax[0].plot(wl, spec2, "g", label="spec2")
    ax[0].legend()
    ax[0].set_ylabel("spectra")
    
    resid1 = spec0 - spec1
    if abs_err:
        resid1 = np.abs(resid1/spec0)
    ax[1].plot(wl, resid1 * 100, "k", label="spec0 - spec1")
    if np.any(spec2):
        resid2 = spec0 - spec2
        if abs_err:
            resid2 = np.abs(resid2/spec0)
        ax[1].plot(wl, resid2 * 100, "b", label="spec0 - spec2")
    ax[1].legend()
    ax[1].set_xlabel(r"$\lambda$\AA")
    ax[1].set_ylabel("residuals (\%)")

    plt.show()

In [8]:
def create_resid_spec(spec0, spec1):
    '''
    Given a base spectrum and a spectrum slightly offset, return the residual spectrum between the two
    '''
    #Normalize both spectra
    spec0 = spec0/np.mean(spec0)
    spec1 = spec1/np.mean(spec1)
       
    resid = np.abs((spec0 - spec1)/spec0)
    return resid

In [9]:
def create_min_spec(base, temp_spec, logg_spec, Z_spec):
    '''
    Given a base spectrum and three perturbing spectra, determine for each pixel the minimum offset so that we can use 
    this as an error envelope in the approximate spectra.
    '''
    #Normalize all three spectra
    base = base/np.mean(base)
    temp_spec = temp_spec/np.mean(temp_spec)
    logg_spec = logg_spec/np.mean(logg_spec)
    Z_spec = temp_spec/np.mean(Z_spec)
    
    #Calculate the residual spectrum for each
    resid_temp = np.abs((base - temp_spec)/base)
    resid_logg = np.abs((base - logg_spec)/base)
    resid_Z = np.abs((base - Z_spec)/base)
    
    #For each pixel, take the smallest value
    #Vstack the arrays and take the min along axis=1
    arr = np.vstack((resid_temp, resid_logg, resid_Z))
    resid_min = np.min(arr, axis=0)
    
    return resid_min

In [10]:
myModelHA.update_Model({"temp":4000, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.})
base = myOrderModelHA.get_spectrum()

myModelHA.update_Model({"temp":4010, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.})
temp_spec = myOrderModelHA.get_spectrum()

myModelHA.update_Model({"temp":4000, "logg":4.05, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.})
logg_spec = myOrderModelHA.get_spectrum()

myModelHA.update_Model({"temp":4000, "logg":4.00, "Z":-0.05, "vsini":0., "vz":0, "logOmega":0.})
Z_spec = myOrderModelHA.get_spectrum()

In [12]:
min_spec = create_min_spec(base, temp_spec, logg_spec, Z_spec)

In [13]:
plt.plot(wl, min_spec * 100)
plt.show()

Temperature


In [14]:
params = {"temp":4000, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux0 = myOrderModel.get_spectrum()


interpolated fl is len 5509

In [15]:
acc_spec = create_resid_spec(base, model_flux0)

In [16]:
fig, ax = plt.subplots(nrows=2, figsize=(11,8), sharex=True)

ax[0].plot(wl, base, "b", label="base")
ax[0].plot(wl, model_flux0, "r", label="approx")
ax[0].legend()
    
ax[1].plot(wl, acc_spec, "r", label="approx")
ax[1].plot(wl, min_spec, "k", label="errenvelope")
ax[1].legend()
plt.show()

In [7]:
params = {"temp":6010, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux1 = myOrderModel.get_spectrum()


interpolated fl is len 5509
before vsini
fl len 5509
5509
self.wl is len 5509
[ 5129.98582119  5130.00574617  5130.02567122 ...,  5240.8739723
  5240.89432797  5240.91468371]
[ 10206524.4374997   10211517.12499971  10205949.31249968 ...,
   9740877.31249984   9782042.81249979   9828856.49999984]
5509 5509

In [8]:
params = {"temp":6000, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModelHA.update_Model(params)
model_flux0HA = myOrderModelHA.get_spectrum()

In [31]:
plot_diff_spec(wl, model_flux0, model_flux0HA, abs_err=True)

In [24]:
plot_diff_spec(wl, base, model_flux0, abs_err=True)

In [ ]:

logg


In [124]:
params = {"temp":2300, "logg":4.05, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux2 = myOrderModel.get_spectrum()

In [127]:
plot_diff_spec(wl, model_flux0, model_flux2, abs_err=True)

In [126]:
plot_diff_spec(wl, model_flux0, model_flux1, model_flux2, abs_err=True)

Metallicity


In [128]:
params = {"temp":2300, "logg":4.00, "Z":-0.05, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux3 = myOrderModel.get_spectrum()

In [130]:
plot_diff_spec(wl, model_flux0, model_flux3, abs_err=True)

In [95]:
plot_diff_spec(wl, model_flux0, model_flux2, model_flux3, abs_err=True)

Using TRES spectra

Perturbations about T=6000, logg=4.0, Z = 0.00

Measured in order 23

It seems like at this stage, with the Mg b lines, that metallicity contributes to the greatest change in the spectrum.

  • Temperature is 0.7%
  • logg is 1.8%
  • Metallicity is 3.7%

So this is interesting, a 1.5% change in flux level corresponds to a 0.05 dex change in metallicity.

Perturbations about T=4000, logg=4.0, Z= 0.0

Measured in order 23

  • Temperature change is about 2.0%
  • Logg change is 2.5%
  • Metallicity change is 5%

Perturbations about T=2300, logg=4.0, Z= 0.0

Measured in order 23

  • Temperature change is 10%
  • Logg change is 7%
  • Metallicity change is 10%

Measured in order 41

  • Temperature change is 4%
  • Logg change is 5%
  • Metallicity change is 5%

Checking that approximate schemes still work

The previous difference spectra were determined using our most accurate implementation of the spectrum modifier yet, that means using the full-resolution PHOENIX grids and k=5 polynomials, such that the spectra are accurate to within floating point precision. Since we have interpolated near a grid point, this also means that any linear approximation to the spline interpolation should be accurate enough.

To test any new approximate scheme for generating spectra, we want to see whether the residual spectrum will be bounded by the smallest spectrum in the {temp, logg, Z} trio, meaning that any error in the approximation will at most contribute to 10% error in any of {temp, logg, Z}. We can increase this threshold if it makes sense.


In [ ]: